﻿#region Copyright information
// 
// Copyright © 2005-2013 Yongkee Cho. All rights reserved.
// 
// This code is a part of the Biological Object Library and governed under the terms of the
// GNU Lesser General  Public License (LGPL) version 2.1 which accompanies this distribution.
// For more information on the LGPL, please visit http://bol.codeplex.com/license.
// 
// - Filename: FDistribution.cs
// - Author: Yongkee Cho
// - Email: yongkeecho@gmail.com
// - Date Created: 2012-09-06 11:39 AM
// - Last Modified: 2013-01-25 3:59 PM
// 
#endregion
using System;
using System.Collections.Generic;
using BOL.Maths.Functions;

namespace BOL.Maths.Distributions
{
    /// <summary>
    /// Represents a F distribution with d1 and d2, degrees of freedom.
    /// <see>Mathematica: ref/FDistribution</see>
    /// <see>http://en.wikipedia.org/wiki/F-distribution</see>
    /// </summary>
    public sealed class FDistribution : IUnivariateDistribution<double>, IEquatable<FDistribution>
    {
        #region Private variables

        private int _d1, _d2;

        #endregion

        #region Public properties

        /// <summary>Gets or sets the number of degrees of freedom, d1 of a F distribution.</summary>
        public int D1 { get { return _d1; } set { _d1 = value; } }

        /// <summary>Gets or sets the number of degrees of freedom, d2 of a F distribution.</summary>
        public int D2 { get { return _d2; } set { _d2 = value; } }

        /// <summary>Gets the range of a F distribution.</summary>
        public IRange<double> Domain { get { return new Range<double>(0, double.PositiveInfinity); } }

        /// <summary>Gets the mean of a F distribution.</summary>
        public double Mean
        {
            get
            {
                if (_d2 < 2)
                    throw new Exception("Mean is defined for d2 >= 2.");

                return _d2/(double) (_d2 - 2);
            }
        }

        /// <summary>Gets the median of a F distribution.</summary>
        public double Median { get { throw new Exception("Median is undefined for F distribution."); } }

        /// <summary>Gets the mode of a F distribution.</summary>
        public double Mode
        {
            get
            {
                if (_d1 < 2)
                    throw new Exception("Mean is defined for d2 >= 2.");

                return (_d1 - 2) * _d2 / (double)_d1 / (_d2 + 2);
            }
        }

        /// <summary>Gets the variance of a F distribution.</summary>
        public double Variance
        {
            get
            {
                if (_d2 < 4)
                    throw new Exception("Mean is defined for d2 >= 4.");

                return 2 * _d2 * _d2 * (_d1 + _d2 - 2) / (_d1 * Math.Pow(_d2 - 2, 2) * (_d2 - 4));
            }
        }

        /// <summary>Gets the skewness of a F distribution.</summary>
        public double Skewness
        {
            get
            {
                if (_d2 < 6)
                    throw new Exception("Mean is defined for d2 >= 6.");

                return (2 * _d1 + _d2 - 2) * Math.Sqrt(8 * (_d2 - 4)) / (_d2 - 6) / Math.Sqrt(_d1 * (_d1 + _d2 - 2));
            }
        }

        /// <summary>Gets the kurtosis of a F distribution.</summary>
        public double Kurtosis { get { throw new Exception("Kurtosis is undefined for F distribution."); } }

        /// <summary>Gets the entropy of a F distribution.</summary>
        public double Entropy { get { throw new Exception("Entropy is undefined for F distribution."); } }

        #endregion

        #region Constructor

        /// <summary>
        /// Instantiates a F distribution.
        /// </summary>
        /// <param name="d1">degree of freedom 1</param>
        /// <param name="d2">degree of freedom 2</param>
        public FDistribution(int d1, int d2)
        {
            if (d1 <= 0)
                throw new ArgumentOutOfRangeException("d1", "d1 is expected to be a positive integer.");
            if (d2 <= 0)
                throw new ArgumentOutOfRangeException("d2", "d2 is expected to be a positive integer.");
            
            _d1 = d1;
            _d2 = d2;
        }

        #endregion

        #region ICloneable implementation

        public FDistribution Clone()
        {
            return new FDistribution(_d1, _d2);
        }

        object ICloneable.Clone()
        {
            return Clone();
        }
        
        #endregion

        #region Public methods

        /// <summary>
        /// Returns the F density given value.
        /// </summary>
        /// <param name="value"></param>
        /// <returns>F density</returns>
        public double Pdf(double value)
        {
            var numerator = Math.Sqrt(Math.Pow(_d1 * value, _d1) * Math.Pow(_d2, _d2) / Math.Pow(_d1 * value + _d2, _d1 + _d2));
            var denominator = value * DistributionFunctions.Beta(_d1 / 2.0, _d2 / 2.0);
            
            return value > 0 ? numerator/denominator : 0;
        }

        /// <summary>
        /// Returns the cummulative F density given value.
        /// </summary>
        /// <param name="value"></param>
        /// <returns>cummulative F density</returns>
        public double Cdf(double value)
        {
            return DistributionFunctions.IncompleteBeta(_d1 / 2.0, _d2 / 2.0, _d1 * value / (_d1 * value + _d2));                                    
        }

        /// <summary>
        /// Returns the value of the F distribution for given probability.
        /// </summary>
        /// <param name="p">probability</param>
        /// <returns></returns>
        public double Quantile(double p)
        {
            if (p.Equals(0))
                return 0.0;

            if (p.Equals(1))
                return Double.MaxValue;

            var beta = new BetaDistribution(_d1 * .5, _d2 * .5).Quantile(p);
            return beta < 2.23e-308 ? Double.MaxValue : _d2 * beta / (1.0 - beta) / _d1;
        }

        /// <summary>
        /// Returns a random value from the F distribution.
        /// </summary>
        /// <returns>random value</returns>
        public double Sample(Random r)
        {
            throw new NotImplementedException();
        }

        /// <summary>
        /// Estimates parameters of the F distribution using maximum likelihood.
        /// </summary>
        public void MaximumLikelihoodEstimate(IEnumerable<double> source)
        {
            throw new NotImplementedException();
        }

        /// <summary>
        /// Estimates parameters of the F distribution using Bayesian methods.
        /// </summary>
        public void BayesianEstimate(IEnumerable<double> source, params dynamic[] priorParameters)
        {
            throw new NotImplementedException();
        }

        #endregion

        #region IEquatable<FDistribution> implementation

        public bool Equals(FDistribution other)
        {
            return _d1.Equals(other._d1) & _d2.Equals(other._d2);
        }

        #endregion

        #region Object overriden

        public override int GetHashCode()
        {
            return D1.GetHashCode() ^ D2.GetHashCode();
        }

        public override bool Equals(object other)
        {
            if (other == null)
                throw new ArgumentNullException("other");

            if (!(other is FDistribution))
                throw new InvalidCastException("The 'other' argument is not a FDistribution object.");
            
            return Equals(other as FDistribution);
        }

        public override string ToString()
        {
            return String.Format("F(d1 = {0}, d2 = {1})", _d1, _d2);
        }

        #endregion
    }
}
